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Abstract. We numerically study the collective coherent and dissipative dynamics in spin lattices with long 
range interactions in one, two and three dimensions. For generic geometric configurations with a small 
spin number, which are fully solvable numerically, we show that a dynamical mean-field approach based 
upon a spatial factorization of the density operator often gives a surprisingly accurate representation of 
the collective dynamics. Including all pair correlations at any distance in the spirit of a second order 
cumulant expansion improves the numerical accuracy by at least one order of magnitude. We then apply 
this truncated expansion method to simulate large numbers of spins from about ten in the case of the 
full quantum model, a few thousand, if all pair correlations are included , up to several ten-thousands 
in the mean-field approximation. We find collective modifications of the spin dynamics in surprisingly 
large system sizes. In 3D, the mutual interaction strength does not converge to a desired accuracy within 
the maximum system sizes we can currently implement. Extensive numerical tests help in identifying 
interaction strengths and geometric configurations where our approximations perform well and allow us to 
state fairly simple error estimates. By simulating systems of increasing size we show that in one and two 
dimensions we can include as many spins as needed to capture the properties of infinite size systems with 
high accuracy. As a practical application our approach is well suited to provide error estimates for atomic 
clock setups or super radiant lasers using magic wavelength optical lattices. 

PACS. 42.50.-p Quantum Optics 06.20.-fMetrology 37.10.JkAtoms in optical lattices 


1 Introduction 

Ensembles of interacting spins in various geometries have 
been at the heart of quantum statistical physics since the 
first models on magnetism were proposed [T] . As the spin- 
spin interaction is nonlinear and the corresponding Hilbert 
space grows exponentially with the number of spins, exact 
analytic as well as full numeric solutions are only possi¬ 
ble for very special cases and geometries [ll[31ll] or a small 
number of spins. The complexity increases even further 
for an open system including collective decay. As a very 
successful approximate numerical approach based on the 
factorization of single site expectation values, a dynam¬ 
ical mean-field method was developed for efficient treat¬ 
ment of larger systems [S] . On the one hand it allowed for 
analytical results in the large dimensions limit while, 
on the other hand, it soon proved very useful for a nu¬ 
merical treatment in low dimensions. Subsequently, the 
general idea of the method was successfully applied to a 
wide range of solid state physics models in the very low 
temperature quantum domain [7] . Recently, this approach 
has proven useful in the description of ultra-cold particle 
dynamics in optical lattices mi- 

The present work is motivated by another, more re¬ 
cent implementation of spin lattices based on ultra-cold 


atoms or molecules trapped in an optical lattice, which 
nowadays can be prepared almost routinely in the labo¬ 
ratory with well defined filling factors and close to zero 
temperature mm- When excited on an optical or in¬ 
frared transition, the trapped particles will interact via 
dipole-dipole energy exchange forming collective excita¬ 
tions mm- In addition, optical transitions intrinsically 
exhibit dissipation via spontaneous decay, which in such a 
lattice becomes a collective effect leading to super- or sub¬ 
radiance [13]. In order to consistently treat such an open 
system, one has to start from a master equation instead 
of the Schroedinger equation after having traced out the 
electromagnetic vacuum modes m- While the interaction 
between a pair of spins can be rather small at a larger dis¬ 
tance, the collective effect of a sizable number of particles 
can still generate noticeable effects in this case|13j. 

Besides using polar molecules, which can possess rela¬ 
tively strong dipole moments m, another interesting im¬ 
plementation is based on using long lived atomic clock 
transitions in a differential light-shift-free magic wave¬ 
length lattice HZldH], where one obtains extremely well 
controllable and precisely measurable systems to study 
even weak spin interactions m and collective decay via 
dipole-dipole energy exchange M- For sufficient densities 
the particles’ effective transition frequency and sponta- 
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neous decay is modified by dipole-dipole interaction [ 20 ] . 
which in turn will influence the performance of a corre¬ 
sponding clock or super-radiant laser pT] , 

While the extremely small dipole moment of a clock 
transition keeps these interactions small, even tiny shifts 
and broadenings will ultimately influence clock accuracy 
and precision. Hence, reliable and converging numerical 
models are required to estimate these effects to many dig¬ 
its, in particular as one tries to work with as large as pos¬ 
sible an ensemble to reduce measurement time and pro¬ 
jection noise. For rather small atom numbers, up to about 
10 , a numerical solution of the full master equation is still 
possible PO] and has showed that shifts and broadening 
can be non-negligible. For larger ensembles at low densi¬ 
ties a so called cluster approach based on statistical aver¬ 
aging of important small particle number configurations 
has already produced first estimates of their scaling with 
the system’s size [22] • Here, we focus on a more generally 
valid approach, namely the above mentioned mean-field 
plus pair-correlation method (MFC) to tackle large sys¬ 
tems at high densities, i.e. up to unit filling. The long range 
nature of the dipole coupling is accounted for by adding 
higher order corrections to the standard local factoriza¬ 
tion approach. In particular, for cavity mediated dipole 
interactions or coupling via nano-fibers even infinite range 
interactions have to be considered |23|. The focus of this 
work is put on developing the appropriate general numer¬ 
ical framework to treat such extended open spin lattices 
in various configurations and test their accuracy and con¬ 
vergence properties by means of the example of collec¬ 
tive decay of a highly excited spin state. This should be 
the basis of future more specific work on concrete imple¬ 
mentations of lattice clock Ramsey spectroscopy Plj and 
super-radiant laser setups mM- 

This work is organized as follows. First, we give a short 
description of the system of coupled spins and introduce 
the corresponding master equation governing their time 
evolution including decay. Using generalized factorization 
assumptions for the density operator of the system we de¬ 
rive two approximate and numerically advantageous meth¬ 
ods to calculate its time evolution. In the subsequent sec¬ 
tion we perform calculations for large ensembles studying 
the influence of long-range interactions. Here, we add an 
extensive numerical analysis to characterize the magni¬ 
tude and scaling of the error of these two approximations 
depending on the geometry and the choice of initial state. 
Finally, we use this method to simulate systems of increas¬ 
ing size to study to which extend a finite sized sample can 
capture the dynamics of larger or even infinite systems. 

2 Interacting spin dynamics 

We consider a system consisting of N two-level subsys¬ 
tems with transition frequency ojq and decay rate 7 in an 
arbitrary spatial configuration. Each particle couples to 
the modes of the free electro-magnetic field and therefore 
all particles are indirectly coupled to one another. Math¬ 
ematically, this problem can be simplified by treating the 
electro-magnetic modes as a single bath and introducing 


effective particle-particle interactions and effective decay 
of particle excitations into this bath, according to m- 
The time evolution of the N spins in a rotating frame 
corresponding to is then governed by a master 

equation 

P =+C[p] (1) 

with the Hamiltonian 

H = ( 2 ) 

i3-,i¥=3 

and Lindblad-term 

The dipole-dipole interaction fiij = ^^G{kQrij) and the 
collective decay Fij = ^jF{kQrij) can be obtained analyt¬ 
ically with 

sin^ ^/cos^ sinA , 

m = + /3 (4a) 

cos^ ^/sin^ cosA , 

G{0 = + /3 (4b) 

with a = 1 —cos^ 9 and (3 = 1 —3 cos^ 6, where 6 represents 
the angle between the line connecting atoms i and j and 
the common atomic dipole orientation m- 

While in systems consisting of only very few particles 
we can study the time evolution by directly integrating the 
master equation, the exponential scaling of the dimension 
of the Hilbert space soon defeats any numerical abilities. 
To be able to represent the state of such a high particle 
number system in a computer one has to make simpli¬ 
fying assumptions about the form of the density matrix. 
In our following calculations we will truncate correlations 
between the particles at a certain order which greatly re¬ 
duces the space needed to store the state of the system in 
memory and allows for treating larger particle numbers. 


2.1 Mean-field method: product state assumption 

In the first nontrivial approximation we neglect correla¬ 
tions altogether and assume that the system is at all times 
in a product state of the subsystems at each site. The den¬ 
sity matrix is approximated by p = pG ), which is also 

called mean-field approximation. The time evolution of the 
system is then governed by the local on site density ma¬ 
trices, which for two-level systems can be obtained from 
a complete set of expectation values for each spin, i.e. the 
expectation values of the Pauli operators {(Jx), (cry) and 
(oz) for a spin 1/2 system. Using this Pauli representa¬ 
tion we need three real numbers to characterize the state 
of each of the two-level sub-systems at a certain point in 
time. The resulting equations for the local spin provide an 
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intuitive insight into the corresponding physics. Explicitly 
we get: 

K) = E ^ E 

i;i^k i;i^k 

(5a) 

(o-fe) = - E ^ E 

i\i^k i;i^k 

(5b) 

(al) = -z E +7(1 - K)) 

i\i^k 

+ ^ E + {a^al)'j (5c) 

i;i^k 

These equations still contain two-particle expectation val¬ 
ues of the form {crfa^), which according to our above as¬ 
sumption can be factorized, i.e. ~ (cf )(o'^). As we 

will see in the next section for weak inter-particle interac¬ 
tions this gives a surprisingly good approximation to the 
interaction induced shifts and can also account for spatial 
inhomogeneities of the system. 


2.2 Extended mean-field method including 
pair-correlations (MPC) 

As a next-order correction to the above mean-field ap¬ 
proach we now include pair-correlations but still neglect 
all higher-order correlations. To this end the density ma¬ 
trix can be approximated by p = 0 ^ pb) ^pb+) 0 

where the first term is the previously used 

product state and the correlations are captured in the op¬ 
erators pb+). The correlations thus have to be chosen 
to generate vanishing single particle expectation values, 
i.e. Tr|(T“pb+)| = 0. Deriving the equations of motion 
in terms of expectation values of Pauli operators leads to 
the same equations as in the mean-field case (eq. . The 
two-particle expectation values are then determined via 
a set of additional equations for the expectation values of 
all two-particle Pauli operator pairs of the type (crfcr^). In 
principle, there are nine such quantities for any pair of par¬ 
ticles pb+). For symmetry reasons three of them are triv¬ 
ially obtained from the others. Similarly to the mean-field 
in the equations for these two-particle correlations higher 
order three-particle correlations appear, which based on 
our assumption of the form of the density operator are 
again approximated by 

{aylal) « -2{at){al){al) + (a“)(a>,") + (af)(a“a^) 

+ ( 6 ) 

Although the resulting equations of motions for the two- 
particle correlations are arguably bulky, we want to dis¬ 
play them explicitly, as the form an essential basis of our 


work. 

{^k^i)= E + E 

- j{alaf) + - i(af)) 

-I E E 

j\j¥^k,l j\j9^k,l 

(7a) 

K^f) = - E E 

r,j¥^k,i j;j¥=k,i 

-7(a|af) + A,((a^af) - i(aD - ^(af)) 

E ^kjifTkafa]) E rijialafa^) 

r,3¥^k,i r,j¥^k,i 

(7b) 

) = E )) 

j\j¥^k,l 

j-,j¥=k,l 

-27(a,X)+7((^f) + (^D) 

+ A/((o-fcCrf) + (c^feCrf)) 

+ 2 E rk,{{alata^) + {ayafay)) 

r,i¥^k,i 

+ ^ E ri,[{alata^) + {alayay)) (7c) 

j'd¥^k,i 

- (erf)) + E ^kik(yla\afi 

r,3=fk,l 

- E f?y(rTfafaf)-7(afaf) 

r,j¥^k,i 

E E ri.ialata)) 

j\i¥^k,l j-,ji^k,l 

(7d) 

(^fc^f) = + E ^kj{<yl<yi(T]) 

j\j¥=k,l 

+ E f?/,((afafaj)-(afara})) 

j-,j¥^k,l 

- - Az((rTfar) - \{at)) 

-I E 

r,3=fk,l 

+ l E ^h-((^k^r^3} + (^k^!^!}) (7e) 

r,3^k,l 

{aiat) = -nu{at)- E ^k,{olalo^) 

3\3=fk,l 

+ E 

3\3¥^k,l 
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- + l{^l) - Ai(KCTf) - 

r,j¥=k,i 





Note, that the number of equations to be solved increases 
quadratically with the number of particles, as we include 
all possible two-particle combinations. This is exponen¬ 
tially slower than the growth of the corresponding Hilbert 
space. In many cases one might even be able to restrict 
this to nearest neighbor couplings only, but for long range 
dipole or cavity mediated interactions, in which we are in¬ 
terested here, no such truncations can be performed safely. 
In principle the method, which in many respects resem¬ 
bles the known cumulant expansion method |26j . can be 
extended towards higher order. However, as we will see 
below, it is already very accurate for our purposes so that 
we will not pursue this task. 

3 Numerical accuracy of the mean-field 
method and second order corrections 

In the previous section we have presented two numerical 
approaches to approximate the master equation (eq. ^ 
by neglecting higher-order quantum correlations. In order 
to examine for which conditions these assumptions lead 
to accurate solutions, we compare these approximations 
with the numerical solution of the full master equation 
for different spatial arrangements, numbers of particles 
and initial states. Additionally we also calculate the case 
of independent particles, which allows us to identify exam¬ 
ples where the error of the approximation is small due to 
a negligible influence of the dipole-dipole-interaction and 
the collective decay. 


3.0.1 Spin dynamics 

To obtain a first intuitive understanding for the quality 
of the different methods we compare the time evolution 
of the expectation values of the Pauli operators for three 
different geometries, i.e. a chain (fig. [^, a square lattice 
(fig. [|) and a cube (fig. [^. As a generic physical exam¬ 
ple, we start with a product state of all spins pointing in 
x-direction. This is the state prepared in the first step of 
a typical Ramsey spectroscopy procedure. It is fully su¬ 
perradiant when all particles are confined in a very small 
spatial volume. Clearly, the dynamics of all three cases 
is significantly different, but they all share certain fea¬ 
tures. First, the solution of the full master equation devi¬ 
ates drastically from the independent particle case, which 
means that the effect of the collective interaction is signif¬ 
icant. This deviation is almost perfectly captured by the 
second order MPC solution, which is, at least visually, al¬ 
most identical to the full solution of the master equation. 


Surprisingly, the mean-field solution shows a qualitatively 
similar behavior already, although it is noticeably not as 
accurate. Note, that for the case of the cube (fig. both 
methods predict the subradiance of the initial spin state 
well |24j . Let us now turn from a visual to a more system¬ 
atic numerical error estimation. 


3.0.2 Systematic accuracy analysis 

In the following we will perform a more rigorous, quanti¬ 
tative analysis for a large range of parameters. In order to 
do this effectively we need a simple measure of accuracy of 
the different methods. A frequently used tool, especially 
in quantum information, is the trace distance which is de¬ 
fined as T(p, cr) = i|Ai| where the Xi are the eigenvalues of 
the matrix representation of p—a. For qubits this measure 
has a very intuitive interpretation, it is just half of the ge¬ 
ometric distance of the two states on the Bloch sphere. In 
fig. I^we use this trace distance between the solution of the 
master equation and the previously presented numerical 
methods at equal points in time to characterize the error 
of the different approximations. In all our examples we ini¬ 
tially start in a product state, which means that the error 
at t = 0 is always zero and since no additional pumping 
is included the system decays to the ground state and the 
trace distance in the long time limit will vanish for all nu¬ 
merical methods. Instead of inspecting the variation of the 
trace distance over time we will use the time-maximum of 
the trace distance as a characterization of the error. 


3.1 Geometry dependence 

In this section we study the geometry dependence of the 
error of the numerical methods measured by the previ¬ 
ously introduced time maximum of the trace distance. We 
distinguish between systems of different dimensionality, a 
ID chain consisting of 8 particles (fig.[^, a 3x3 section of a 
2D square lattice (fig.[^ and a cube as a 3D configuration 
(fig-0). For each of these examples we calculate the depen¬ 
dence of the error on the distance between the particles. 
Further on we vary the initial state and the orientation of 
the polarization vector and show three typical results. In 
the following the initial state is characterized by O which 
measures the polar angle between a given state towards 
the ground state on the Bloch sphere. Several interesting 
features stand out immediately. The bigger the distance 
between the particles, the smaller the error of neglecting 
higher-order correlations. As can be seen from the trace 
distance between the solution of the master equation and 
the independently decaying case this is to some degree an 
artifact of the decreasing strength of the dipole-dipole in¬ 
teraction which in the far field limit has a ^ dependency 
but at least for MPC the error decreases much faster. In 
nearly all cases the mean-field approach yields a noticeable 
improvement, yet, when all spins start in the excited state 
it reproduces the results of independent particles only. In 
fact, as one can show from the mean-field equations, in this 
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Fig. 1. Time evolution of the expectation values of the Pauli operators Oy and of the central spin in a chain consisting 
of 7 spins with spin-spin distance d — O.ISAq. The system is simulated using independent spins (red), the mean-field method 
(blue), MPC (green) and by solving the whole master equation (dashed black). The dipole is orientated orthogonally to the 
chain. 
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Fig. 2. Time evolution of the expectation values of the Pauli operators CTj,, Oy and Oz of the central spin in a 2D-square lattice 
consisting of 3x3 spins with nearest spin-spin distance d = 0.5Ao. The system is simulated using independent spins (red), the 
mean-field method (blue), MPC (green) and by solving the whole master equation (dashed black). The dipole is orientated 
orthogonally to the plane. 





Time [7 Time [7 Time [7 

Fig. 3. Time evolution of the expectation values of the Pauli operators CTj,, (jy and Oz of a single spin in a cube configuration 
with nearest spin-spin distance d — O.6A0. The system is simulated using independent spins (red), the mean-field method (blue), 
MPC (green) and by solving the whole master equation (dashed black). The dipole is orientated orthogonally to an arbitrary 
face of the cube. 


case the time evolution is completely identical, so mean- 
field results in no improvement over simply ignoring the 
collective effects. 

3.2 Initial state dependence 

To further analyze the dependence of the error on the ini¬ 
tial state we consider a chain of six particles with three 
different particle distances. Initially the system is in a 
product state where all single particles are in the same 
Bloch state. For simplicity we only consider pure states 


and since the time evolution is invariant under a global 
rotation around the z-axis the only remaining variable is 
the polar angle 0. In fig. the dependence of the error 
on this polar angle is shown. For 6> = 0 the system is in 
the ground state and the error vanishes. For a small exci¬ 
tation the mean-field method gives a substantial improve¬ 
ment compared to the independently decaying system but 
for a nearly totally excited state the advantage disappears 
more and more. In contrast, MPC performs convincingly 
for all initial states. 
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Time [7 Time [7 Time [7 

Fig. 4. Trace distance between the density operators calculated by the master equation and density operators calculated for 
independent spins (red), mean-field (blue) and MFC (green) for the previously discussed chain configuration (a), 3x3 square 
lattice (b) and cube configuration (c). 




d/Ao d/Ao d/Ao 

Fig. 5. Distance dependency of the time-maximum of the trace distance between the results of the master equation and the 
results of an independent evolution (red), mean-field (blue) and MFC (green) for a chain consisting of 8 spins for different initial 
states and dipole orientations, (a) O = 7r/2, Cdipoie = e^. (b) 6> = tt, Cdipoie = Cz- (c) 0 = •7r/2, Cdipoie = e^,. 



Fig. 6. Distance dependency of the time-maximum of the trace distance between results of the master equation and results of 
independent evolution (red), mean-field (blue) and MFC (green) for a square lattice consisting of 3x3 spins for different initial 
states and dipole orientations, (a) O = 7r/2, Cdipoie = e^. (b) 0 = tt, Cdipoie = e^. (c) 0 = •7r/2, Cdipoie = e^,. 


3.3 Spin-number dependence 

Finally, we investigate the dependence of the error on the 
number of particles in the system, i.e. a chain consist¬ 
ing of N particles. The result of this analysis is shown in 
fig. |§i. In this double logarithmic plot the error appears 
to be nearly linear but slightly shifted for varying particle 
numbers which leads us to the following estimate for the 
error 

err{N,d) = CN*d^^. ( 8 ) 

The exponent fcjv and the factor Cn can be extracted 
from this error plot and are shown in fig. and fig. 


respectively. The error exponent turns out to be indepen¬ 
dent of the number of particles and is -1 for independently 
decaying spins which is not surprising since the collective 
interaction in the far field drops with i. However, increas¬ 
ing the distances doesn’t improve the mean-field results 
whereas MFC has an error exponent of -2 and gains dras¬ 
tically on accuracy. 

4 Approximation of very large (infinite) 
systems 

Recent research on the effect of geometry on the per¬ 
turbation of the spin dynamics by collective interactions 


































































S. Kramer and H. Kitsch: Generalized mean-field approach 


7 



d/Ao 




d/Ao d/Ao 


Fig. 7 . Distance dependency of the time-maximum of the trace distance between the results of the master equation and the 
results of independent evolution (red), mean-field (blue) and MFC (green) for 8 spins in a cube configuration for different initial 
states and dipole orientations, (a) 0 = 7r/2, Cdipoie = e^. (b) 6 > = tt, Cdipoie = Cz- (c) 0 = •7r/2, Cdipoie = + Gy-\- Cz). 





Fig. 8. Dependence of the time-maximum of the trace distance between the results of the master equation and the results of 
an independent evolution (red), mean-field (blue) and MPC(green) on the initial Bloch state characterized by the polar angle 
0 for a chain consisting of 6 spins with spin-spin distance d = 0.5Ao (a), d — l.OAo(b) and d = lO.OAo (c). 





Fig. 9. Dependence analysis of the time-maximum trace distance between results of the master equation and results of an 
independent evolution (red), mean-field (blue) and MFC (green) on the spin distance of a chain consisting of = 3... 9 spins. 
Higher spin numbers correspond to slightly increased trace distances (a). An approximation of the trace distances by Cn * 
results in the spin-number dependency of the error-exponent fcjv (b) and the error-factor Cn (c). 


was mostly limited to systems consisting of only very few 
atoms. In lack of better alternatives one might be tempted 
to extrapolate results obtained from these small-sized sys¬ 
tems to larger ensembles but in general this attempt could 
fail miserably. Armed with the knowledge about the accu¬ 
racy of the mean-field and MFC methods and their ability 
to simulate moderately large systems we can use them to 
investigate how many particles are needed to make sat¬ 
isfying statements about infinite systems. More precisely, 
we want to know how collective spin quantities of the type 
will change for different numbers of particles. 
Of course, it is not a priori clear if these expectation val¬ 


ues will converge at all. To answer this question we will 
study two different examples. 


4.1 Linear equidistant chain 


We consider a Wparticle spin chain with particle distance 
d and calculate the dynamics of the whole system where 
initially all spins are in the {ax} = 1 state. Tracing out all 
but the innermost spin allows us to compare the dynamics 
of this single spin for a varying number of surrounding par¬ 
ticles. The result of this analysis for a certain distance d 
after an integration time of 2"f~^ is shown in fig. 10 For- 
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Fig. 10. Expectation values of the Pauli operators {ox), {oy) and (ctz) of the central spin in a chain consisting of N spins with 
distance d — O.QAq after a time evolution for 27 “^. Initially all spins are in the state (a^) = 1 and the system is solved for 
independent spins (red), using mean-field (blue), MPC (green) and the master-equation (dashed black). 



Fig. 11. Time-maximum of the trace-distance between the reduced density matrix of the central spin state in a chain consisting 
of N spins compared to the chain consisting of = 20001 (blue) where both quantities are results of mean-field 

simulations. Approximation of this trace-distance by Cd * (dashed yellow) for different spin-spin distances d = 0.7Ao(a), 
d = l.OAo(b) and d = 3.7Ao(c) 





Fig. 12. Time-maximum of the trace-distance between the reduced density matrix of the central spin in a chain consisting 
of N spins compared to the chain consisting of = 401 (blue) where both quantities are results of MPC simulations for 

different spin-spin distances d — 0.7Ao(a), d = l.OAo(b) and d = 3.7Ao(c). Approximation of this trace-distance by Cd * 
(dashed yellow). 


tunately, all methods yield more or less the same result 
and differ significantly from the independently decaying 
case, indicating that the variation for small particle num¬ 
bers and the ultimate convergence for large systems is not 
a numerical artifact. This result hints at the fact that a 
suitable number of particles will indeed give a usable ap¬ 
proximation of big systems. To consolidate this claim we 
perform a more extensive and quantitative test. What we 
actually would like to test is how much the time evolution 
of the single central spin in a chain consisting of N parti¬ 
cles differs from the time evolution of a spin in an infinite 
chain. However, we are not aware of a method to solve 


the infinite chain exactly which leaves us with the option 
to compare the central spin of a particle chain with 
a chain containing as many spins as numerically possible 
only. In fig. [TT] and fig.[^ the dynamics of a 20001 particle 
mean-field simulation and a 401 particle MPC simulation 
are used as the best possible approximation of an infinite 
chain for three different spin-spin distances, respectively. 
In most cases the addition of further spins affects the cen¬ 
tral spin less and less and is approximately linear in this 
double logarithmic plot, i.e. the trace distance between the 
infinite chain and the N-particle chain for a certain dis¬ 
tance can be estimated by T{N,mi) = C^N^. By fitting 
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d/Ao d/Ao 

Fig. 13. Error exponent kd and Error factor Cd dependence on the spin-spin distance determined in the previous fits Cd * 
using mean-field (blue) and MFC (green) in the case of a chain. 


this function to the numerical results we can determine the 
exponent kd and the factor Cd depending on the distance 
which is plotted in fig. [T^ For nearly all distances both 
the mean-field method as well as the MFC method predict 
that adding further particles has an effect proportional to 
^ only, allowing us to easily estimate the number of parti¬ 
cles needed to approximate the infinite chain dynamics to 
a desired accuracy. However, when the spin-spin distance 
is close to a multiple of the transition wavelength Aq, the 
dynamics of the infinite chain never seems to be captured 
by a finite size approximation. 


4.2 Hexagonal lattice 

With these very encouraging results for a ID chain, let us 
see, if this holds true for higher dimensional geometries 
as well. Unfortunately we failed to obtain convincing re¬ 
sults for 3D cubic lattices, since the number of particles 
needed for convergence of the numerical result turns out to 
very much exceed the possibility of MFC and even of the 
mean-field method. For 2D geometries at least, the mean- 
field method delivers some meaningful, albeit by far not as 
beautiful results. Fig. [T^ shows the numerically obtained 
approximations for the error exponent and the error fac¬ 
tor for the case of a hexagonal lattice, where additional 
particles are added in rings around the central spin. The 
outcome looks rather noisy, probably due to too small a 
choice for the number of particles used as an approxima¬ 
tion of the infinite lattice. It turns out that compared to 
the chain a lot more particles are needed to reliably ap¬ 
proximate an infinite hexagonal lattice, i.e. the influence of 
additional particles reduces the error with approximately 
where this exponent is a rather rough estimate. 

5 Numerical complexity of the different 
methods 

Finally we want to add some considerations on the mem¬ 
ory and CFU requirements of the different methods and 


show that our implementations behave as expected in this 
regard. When solving the master equation the state of 
the system is captured as a density matrix of dimension 
2^^. The time evolution according to a master equation 
is equivalent to a matrix-matrix multiplication and there¬ 
fore has a time complexity of 0{2^^). In the case of mean- 
field a state can be characterized by 3N real numbers and 
according to the mean-field equations (eq. the time 
complexity is then approximately 0{N'^). For the MFC 
method the state consists of one mean-field state and nine 
correlation matrices of the form C“.^ = (crfcr^). Using the 

relation means that we need roughly real 

numbers to represent one MFC state. The time complexity 
is, according to the MFC equations (eq.[^, approximately 
0{N^). The results of this analysis are presented in fig. 


15 



Number of Spins 


Fig. 15. Time needed to integrate a spin chain consisting of 
N spins from 0 to on a single CPU. The solid lines are 
the result of the benchmarks for master (black), MFC (green) 
and mean-field (blue), the dashed lines are the corresponding 
theoretical predictions. 
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Fig. 14. Error exponent kd and error factor Cd dependence on the spin-spin distance determined in the previous fits Cd * 
using mean-field in the case of a hexagonal lattice. 


6 Conclusions and outlook 

We have demonstrated that an effective mean-field method 
with added pair correlations provides a numerically effi¬ 
cient and surprisingly accurate method to simulate open 
spin systems with general non local spin-spin interaction 
and collective decay up to moderately high particle num¬ 
bers and significant interaction strength. Particularizing 
to dipole-dipole interaction and collective spontaneous de¬ 
cay has allowed us to establish a numerical estimate of the 
accuracy and scaling properties of our methods. Further¬ 
more we can show for ID chains that tractable system sizes 
already approach the behavior of infinite systems allowing 
for an estimate of the magnitude of the error due to the 
truncation of the system. For 2D systems the lowest order 
mean-field approach still allows to reach adequate system 
sizes to approximate infinite systems, whereas the scaling 
is unfavorable to accurately approximate infinite 3D sys¬ 
tems. In future work, we plan to apply these methods to 
study collectively enhanced as well as suppressed decay 
in magic wavelength lattices for clock atoms. The simula¬ 
tions should also provide us with predictions of geometries 
and excitation schemes to minimize dipole-dipole induced 
shifts in order to improve the accuracy of atomic clocks. 
Possible approaches would be to analyze different geome¬ 
tries, use initial phase spread rotations and spin squeezing. 
As an interesting extension of this model we also want to 
embed such spin systems inside a cavity and derive cor¬ 
responding mean-field and MPC equations for the arising 
infinite range interactions. This should give us a basis to 
simulate super-radiant lasers for larger ensembles includ¬ 
ing their interaction. Note, that as we are simulating an 
open system anyway, including a finite bath temperature 
will hardly change the complexity of these calculations and 
could be used to identify temperature dependent phase 
transitions in the system. 

This work has been supported by the Austrian Science Fund 
FWF through the SFB F4013 FoQuS. 
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